#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <pthread.h>
#include <papi.h>
#include <emmintrin.h>
#include <sys/timerfd.h>

#ifdef TCOND
#include <unistd.h>
#endif

int events[3] = {PAPI_TOT_CYC, PAPI_L1_TCM, PAPI_L2_TCM};
int eventnum = 3;
long long values[3];
int eventset;
FILE *gnuplotfile;

#ifdef FP
typedef float TYPE;
#elif DP
typedef double TYPE;
#else
typedef int TYPE;
#endif

struct periodic_info
{
	int timer_fd;
	unsigned long long wakeups_missed;
};

static int make_periodic (unsigned int period, struct periodic_info *info)
{
	int ret;
	unsigned int ns;
	unsigned int sec;
	int fd;
	struct itimerspec itval;

	/* Create the timer */
	fd = timerfd_create (CLOCK_MONOTONIC, 0);
	info->wakeups_missed = 0;
	info->timer_fd = fd;
	if (fd == -1)
		return fd;

	/* Make the timer periodic */
	sec = period/1000000;
	ns = (period - (sec * 1000000)) * 1000;
	itval.it_interval.tv_sec = sec;
	itval.it_interval.tv_nsec = ns;
	itval.it_value.tv_sec = sec;
	itval.it_value.tv_nsec = ns;
	ret = timerfd_settime (fd, 0, &itval, NULL);
	return ret;
}

static void wait_period (struct periodic_info *info)
{
	unsigned long long missed;
	int ret;

	/* Wait for the next timer event. If we have missed any the
	   number is written to "missed" */
	ret = read (info->timer_fd, &missed, sizeof (missed));
	if (ret == -1)
	{
		perror ("read timer");
		return;
	}

	/* "missed" should always be >= 1, but just to be sure, check it is not 0 anyway */
	if (missed > 0)
		info->wakeups_missed += (missed - 1);
}

void Print(TYPE *C, int N, char *Name)
{
  int i, j;
  FILE *OutputFile;
  OutputFile = fopen(Name, "w");
  fprintf(OutputFile, "Matrix C:\n");
  for (i = 0; i < N; i++) {
    fprintf(OutputFile, "     ");
    for (j = 0; j < N; j++) {
#ifdef FP
      fprintf(OutputFile, " %f", C[i*N+j]);
#elif DP
      fprintf(OutputFile, " %e", C[i*N+j]);
#else
      fprintf(OutputFile, " %i", C[i*N+j]);
#endif
    }
    fprintf(OutputFile, "\n");
  }
}

struct ThrdArg
{
  TYPE *A; // First  Matrix
  TYPE *B; // Second Matrix
  TYPE *C; // Result Matrix
  int T;   // Thread Number
  int P;   // Part
  int N;   // Size
#ifdef TCONDBAR
  pthread_barrier_t *barr;
#endif
};

void *PapiKernel (void *Args) {
	struct periodic_info info;
        long long i=0;

        // 100ms interval
	make_periodic (100000, &info);
	while (1)
	{
		i+=100;
		/* Do useful work */
		PAPI_read(eventset, values);
		fprintf(gnuplotfile,"%lld\t%lld\t%lld\n", i, values[1], values[2]);
		wait_period (&info);
	}
	return NULL;

}

void *MatMulKernel (void *Args)
{
  int i, j, k;
  struct ThrdArg *Argument = (struct ThrdArg *) Args;
  TYPE *A=Argument->A;
  TYPE *B=Argument->B;
  TYPE *C=Argument->C;
  int T=Argument->T;
  int P=Argument->P;
  int N=Argument->N;
  printf("<=== Thread %i -> Before Matrix Multiplication ===>\n", T);
  for (i = 0; i < P; i++) {
#ifdef TCONDBAR
    if((!(i+1) % (P/4))) {
      pthread_barrier_wait(Argument->barr);
    }
#endif
    for (j = 0; j < N; j++)
      for (k = 0; k < N; k++) {
        C[i*N+j] += 
                     A[i*N+k]*
                     B[k*N+j];
      }
  }
  printf("<=== Thread %i -> After  Matrix Multiplication ===>\n", T);
  //pthread_exit(0);
}

// M is the total number of threads.
void MatMul (TYPE *A, TYPE *B, TYPE *C, int N, int M)
{
  int i, j, status;
  struct ThrdArg Argument[M];
  pthread_attr_t attr;
#ifdef BASE
  pthread_t Thread;
#else
  pthread_t Thread[M+1];
#endif

#ifdef TCOND
  struct timespec CondTime;
  struct timeval  ContTv;
  pthread_cond_t  Condition;
  pthread_mutex_t CondMutex;
#endif

#ifdef TCONDBAR
  pthread_barrier_t barrier1; 
  pthread_barrierattr_t barrattr;
  pthread_barrier_init(&barrier1, &barrattr, M);
#endif

  pthread_attr_init(&attr);
  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

for(i = 0; i < M; i++) {
    (Argument+i)->A = A+i*N*N/M;
    (Argument+i)->B = B;
    (Argument+i)->C = C+i*N*N/M;
    (Argument+i)->T = i;
    (Argument+i)->P = N/M;
    (Argument+i)->N = N;
#ifdef TCONDBAR
    (Argument+i)->barr = &barrier1;
#endif
}

#ifdef SEQ
for(i = 0; i < M; i++) {
  pthread_create(Thread+i, &attr, MatMulKernel, (void*) (Argument+i));
  pthread_join(Thread[i], (void **)&status);
}
#elif TCOND
  pthread_mutex_init(&CondMutex, NULL);
  pthread_cond_init (&Condition, NULL);
  for(i = 0; i < M; i++) {
  	if(i != 0) {
      pthread_mutex_lock(&CondMutex);
      pthread_cond_timedwait(&Condition,&CondMutex,&CondTime);
      pthread_mutex_unlock(&CondMutex);
	  }
    pthread_create(Thread+i, &attr, MatMulKernel, (void*) (Argument+i));
    clock_gettime(CLOCK_REALTIME,&CondTime);
    CondTime.tv_nsec += TIME*1000;
  }

  for(i = 0; i < M; i++)
    pthread_join(Thread[i], (void **)&status);
#elif BASE
  pthread_create(&Thread, &attr, PapiKernel, NULL);
  MatMulKernel((void *)Argument);
  pthread_cancel(Thread);
#else
  pthread_create(Thread+M, &attr, PapiKernel, NULL);
  for(i = 0; i < M; i++)
    pthread_create(Thread+i, &attr, MatMulKernel, (void*) (Argument+i));
  
  for(i = 0; i < M; i++)
    pthread_join(Thread[i], (void **)&status);
  pthread_cancel(Thread[M]);
#endif

}

int main (int argc, char *argv[])
{
  int i, j, k, N, M;
  TYPE *A, *B, *C;

  //papi init code begin
  long long s;
  long long e;
  int retval;
  int eventcode;
  PAPI_option_t opt;

  if(PAPI_VER_CURRENT != PAPI_library_init(PAPI_VER_CURRENT)){
    printf("Can't initiate PAPI library!\n");
    exit(-1);
  }

  if (PAPI_thread_init(pthread_self) != PAPI_OK) {
    printf("Can't thread init!\n");
    exit(-1);
  }

  eventset = PAPI_NULL;
  if(PAPI_create_eventset(&eventset) != PAPI_OK){
    printf("Can't create eventset!\n");
    exit(-3);
  }

  if ( ( retval = PAPI_assign_eventset_component( eventset, 0 ) ) != PAPI_OK ) {
    printf("Can't assign_event_component!\n");
    exit(-3);
  }                

  memset( &opt, 0x0, sizeof ( PAPI_option_t ) );
  opt.inherit.inherit = PAPI_INHERIT_ALL;
  opt.inherit.eventset = eventset;
  if ( ( retval = PAPI_set_opt( PAPI_INHERIT, &opt ) ) != PAPI_OK ) {
    printf("Can't set inherit!\n");
    exit(-3);
  }

  //PAPI_event_name_to_code("ix86arch::LLC_MISSES",&eventcode);
  PAPI_event_name_to_code("MEM_LOAD_MISC_RETIRED:LLC_MISS",&eventcode);
  events[1] = eventcode;
  //PAPI_event_name_to_code("LLC_REFERENCES",&eventcode);
  PAPI_event_name_to_code("CPU_CLK_UNHALTED:THREAD_P",&eventcode);
  events[2] = eventcode;
  for(i=0;i<3;i++) {
    retval = PAPI_add_event(eventset, events[i]);
    if(retval != PAPI_OK){
      printf("error %d\n",retval);
      exit(-4);
    }
  }
  //papi init code end

  if (argc == 3){
    N = atoi(argv[1]);
    M = atoi(argv[2]);
  } else {
    N = 512;
    M = 2;
  }
  A = (TYPE *) calloc(N*N,sizeof(TYPE));
  B = (TYPE *) calloc(N*N,sizeof(TYPE));
  C = (TYPE *) calloc(N*N,sizeof(TYPE));
  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++) {
      A[i*N+j] = (TYPE) rand();
      B[i*N+j] = (TYPE) rand();
    }

  /* create output file in gnuplot format */
  gnuplotfile = fopen("tcm_time.txt", "w");
  fprintf(gnuplotfile,"time(100ms)\tMEM_LOAD_MISC_RETIRED:LLC_MISS\tCPU_CLK_UNHALTED:THREAD_P\n");
  s = PAPI_get_real_usec();
  PAPI_start(eventset);
  MatMul(A,B,C,N,M);
  PAPI_stop(eventset, values);
  e = PAPI_get_real_usec();
  fclose(gnuplotfile);
//  Print(A,N,"mmmthra.txt");
 // Print(B,N,"mmmthrb.txt");
  //Print(C,N,"mmmthrc.txt");

  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++) {
      _mm_clflush((const void *)(A+i*N+j));
      _mm_clflush((const void *)(B+i*N+j));
      _mm_clflush((const void *)(C+i*N+j));
    }

  free(A);
  free(B);
  free(C);

  /*Print out PAPI reading*/
  printf("Wallclock time: %lld usec\n",e-s);
  printf("Total Cycles\t%lld\n", values[0]);
  printf("MEM_LOAD_MISC_RETIRED:LLC_MISS\t%lld\n", values[1]);
  printf("CPU_CLK_UNHALTED:THREAD_P\t%lld\n", values[2]);

  return 0;
}
